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1.  Introduction 


Geman  and  Geman  (1984)  discussed  a  methodology  for  pixel  image  restoration  which 
depended  on  the  idea  of  modelling  the  image  by  a  Markov  random  field.  A  key  feature  of  their 
approach  was  the  possible  placing  of  “edge  elements”  at  “line  sites”  between  pixels  of  the  image. 

In  Geman  and  Geman’s  approach,  a  prior  distribution  for  the  image  is  constructed  by  first 
constructing  a  prior  Gibbs  distribution  for  the  process  of  edge  elements  and  then  specifying  the 
prior  for  the  pixels  themselves  conditional  on  the  edge  process.  In  the  specification  of  the  pixel 
process,  contiguous  pixels  separated  by  a  line  site  at  which  an  edge  element  is  actually  present 
are  not  considered  as  neighbours,  and  so  are  allowed  to  have  quite  different  grey  levels  without 
incurring  any  penalty  in  the  prior  likelihood. 

The  edge  process  idea  corresponds  to  the  notion  that  the  image  is  segmented  into  regions 
over  each  of  which  its  behaviour  is  relatively  homogeneous,  or  at  least  is  not  subject  to  abrupt 
changes;  from  one  region  to  another,  however,  large  differences  in  behaviour  are  possible.  The 
changes  in  behaviour  may  relate  either  to  overall  grey  level  or  to  more  subtle  properties  such 
as  texture. 

In  this  report  we,  shall  focus  attention  on  the  specification  of  the  edge  process,  and  show 
how  various  geometrical  insights  suggest  how  the  prior  Gibbs  distribution  should  be  constructed. 
Our  discussion  will  suggest  relative  costs  for  possible  configurations  somewhat  different  from 
those  proposed  by  Geman  and  Geman  (1984).  In  addition  or$r  scheme  will  provide  methods  for 

dealing  with  rectangular  and  irregular  pixel  patterns.  —  ^ 

/ 

The  present  report  is  unashamedly  speculative  and  theoretical.  Practical  implementation 
and  investigation  of  the  ideas  presented  here  is  in  progress  and  will  be  reported  subsequently. 

2.  The  Gibbs  Log  Likelihood  as  a  Penalty  Function 

The  Gibbs  distribution  approach  constructs  a  prior  likelihood  for  the  edge  process  by  first 
defining  a  set  of  cliques  of  line  sites.  Each  clique  C  consists  of  a  small  set  of  sites;  in  the  Geman 
and  Geman  paper  the  cliques  are  the  collections  of  four  line  sites  with  a  common  vertex.  The 
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Figure  3.1:  Possible  types  of  configuration  for  regular  edge  process. 

The  costs  ascribed  to  these  configurations  by  Geman  and  Geman  (1984)  are  given  in  Table 
3.1. 


Table  3.1:  Costs  of  the  configurations  of  Figure  3.1  used  by  Geman  and  Geman  (1984) 


Type  of  configuration  0  1 
Cost  V  0  2.7 


2  3 

1.8  0.9 


4  5 

1.8  2.7 


We  shall  write  v*  for  the  cost  of  a  configuration  of  type  t,  and  explore  the  consequences  of 


various  choices  of  t;,-. 


1(a)  =  c(0)/c  (^)  =  1  l(a\/2).  To  deal  with  \  <  a  <  1,  define  0O  =  tan-1(2a  -  1)  and  rewrite 

c(8)  =  /i-1V3  sec  do  {cos  #o  cos  9  +  sin  0O  sin#} 

(3.2) 

=  h~  V3  sec  #0  cos(#  -  #0). 

Since  for  |  <  a  <  1  we  have  0  <  #o  <  f ,  it  follows  that,  for  a  in  this  range,  c(8)  has  a  maximum 
at  #0  ami  that  1(a)  =  max  (sec  8q,  sec  -  #o)}.  Hence  1(a)  is  minimized  by  setting  $c=  |. 

The  minimum  value  sec|  =  (4  -  2n/2)^2  «  1.082.  Thus  it  follows  that  the  minimax  score 
1(a)  is  optimized  by  setting  2a  -  1  =  tan  |,  which  implies  that  a  —  ^  (l  +  tan  |)  =  l/v^2. 
If  this  value  of  a  is  used,  then  lines  parallel  to  the  lattice  directions  or  those  at  45°  to  these 
directions  will  cost  the  same  amount  per  unit  length,  while  the  most  expensive  lines  will  be 
those  at  22 to  the  axis  directions,  which  will  cost  about  8%  more.  It  is  interesting  to  note 
that  the  Geman-Geman  value  a  =  2  yields  1(a)  =  2>/2  as  2.83,  a  much  larger  value. 

It  can  also  be  shown,  by  somewhat  tedious  algebra,  that  a  =  l/\/2  also  minimizes  other 
criteria  of  variability  of  c(8),  for  example  the  coefficient  of  variation  of  c(0)  with  9  uniformly 
distributed  over  [0,  J] . 

The  arguments  of  this  section  make  it  possible  to  settle  on  a  charge  for  configurations  of 
types  2  and  3.  Suppose  it  is  intended  to  penalize  boundaries  in  the  underlying  picture  by  an 
amount  (3  per  unit  length.  In  an  ideal  world  we  would  like  to  choose  V2  and  v$  in  (3.1)  to  ensure 
that  c(8)  =  (3  for  all  9.  As  we  have  seen,  this  cannot  be  attained  exactly  for  all  9 ,  but  setting 
v2/v3  =  2-1/2  will  minimize  the  variability  of  c(9)  as  9  varies.  Having  settled  the  ratio  U2/V3, 
it  is  natural  to  choose  vj  to  ensure  that  (27r)-1  c(8)d9  =  (3.  By  simple  algebra,  from  (3.2), 

(2ir)_1  f  c(9)d0  =  4t-1  /  h~lvz  sec  (|)  cos  (9  —  |)  d9 
Jo  Jo 

=  87r-1h-1U3  tan  (|)  =  voh~lk~x 
where  the  constant  k  =  gjr/tan  (5)  ~  0.948. 

It  follows  that  setting  03  =  k/3h  and  v2  =  2~x^kf3h  will  ensure  that,  while  c(9)/(3  is  only 
exactly  1  for  certain  values  of  9 ,  it  will  be  the  case  that  c(8)/(3  lies  between  0.948  and  1.027  for 
all  8  and  furthermore  that  the  average  value  of  c(9)  over  (uniformly  distributed)  9  is  precisely 


Let  tib  be  the  number  of  branches  and  nc  the  number  of  crossings.  It  is  immediate  that 


—  Tifc  rig.  (3.3) 

In  order  to  count  the  number  of  boundary  sections,  notice  that  three  boundary  sections 
meet  at  each  branch  and  four  at  each  crossing.  Thus  the  number  of  ends  of  boundary  sections 
in  3n&  +  4nc,  and  since  each  boundary  section  has  two  ends,  we  have 

ne  =  §n4  +  2nc.  (3.4) 

Substituting  (3.3)  and  (3.4)  into  (3.5)  yields 

nI  =  1  +  \nb  +  »c  (3.5) 

Formula  (3.5)  gives  a  natural  price  to  be  charged  for  branches  and  crossings.  If  it  is  desired  to 
penalize  an  amount  p  for  each  region  in  the  pattern,  then  one  should  charge  \p  for  each  branch 
point  and  p  for  each  crossing.  Of  course,  if  the  edge  configuration  gives  rise  to  regions  that 
are  not  simply  connected,  then  the  number  of  regions  will  no  longer  be  given  by  (3.5),  and  the 
charge  (^nj,  4-  nc)  p  will  have  to  be  considered  in  its  own  right  as  a  penalty  for  the  complexity 
of  the  pattern. 

3.3.  Endings 

A  pattern  made  up  of  disjoint  regions  cannot,  of  course,  have  a  configuration  of  edges 
containing  any  endings  at  all.  Therefore  the  philosophy  that  we  are  adopting  would  naturally 
led  to  an  infinite  charge  for  configurations  of  type  1  in  Figure  1.  However  to  set  any  penalty  value 
to  infinity  leads  to  algorithmic  difficulties  in  using  the  model  in  practice,  because  it  yields  a  prior 
model  for  the  edge  process  under  which  some  configurations  have  probability  zero.  This  violates 
the  condition  of  positive  probability  for  all  configurations  under  which  the  theory  and  practice 
of  Markov  random  fields  is  developed;  see,  for  example,  Geman  and  Geman  (1984,  Section  4). 
In  any  case,  it  seems  excessively  dogmatic  to  exclude  certain  configurations  completely,  since 
there  may  be  good  physical  reasons  for  a  boundary  to  peter  out  in  the  middle  of  a  region. 
Therefore  an  approach  that  is  likely  to  be  more  satisfactory  is  to  ascribe  a  cost  A  to  each  “loose 
end”  in  the  boundary  pattern,  where  A  is  set  to  a  relatively  large  value.  Precisely  how  large 


At  the  edge  of  the  pattern,  there  will  not  be  four  line  sites  meeting  at  each  vertex.  A  great 
advantage  of  our  methodology  is  tht  it  makes  it  possible  very  easily  to  find  the  appropriate 
charge  for  the  possible  configurations  that  arise.  Consider  first  the  case  where  the  window 


p-^|l — ^ —  edlg-e.  ua-a^out 


Figure  3.3:  A  simple  configuration  at  the  edge  of  the  window. 


is  locally  aligned  with  the  pixel  direction,  as  in  Figure  3.3.  The  vertex  A  will  be  associated 
with  a  clique  containing  exactly  one  line  site;  if  the  edge  shown  in  Figure  3.3  is  present  then 
the  point  A  will  be  a  branch,  and  so  the  appropriate  charge  is  \p,  while  if  the  edge  is  absent 
the  appropriate  charge  is  0. 

A  rarer  kind  of  incomplete  clique  arises  as  in  Figure  3.4. 


The  case  where  just  a  single  edge  is  present  will,  as  before,  count  as  a  branch  point  and  so 
should  be  charged  \p,  while  the  configuration  containing  both  edges  contributes  two  boundary 
section  ends  and  so  should  be  charged  p,  provided  the  pixels  numbered  1  and  3  are  not  to  be 
regarded  as  neighbours. 


Although,  in  contrast  with  the  case  of  square  pixels,  there  are  fewer  types  of  configuration 
to  consider,  the  irregularity  of  the  pixels  means  that  it  is  no  longer  necessarily  the  case  that  all 
configurations  of  a  particular  type  should  attract  the  same  penalty. 

The  first  stage  in  the  assignment  of  costs  to  various  configurations  is  to  use  the  same 
arguments  as  in  the  square  lattice  case  to  deal  with  configurations  of  types  0, 1  and  3,  which  are 
charged  0,  A  and  \p  respectively.  It  remains  to  ascribe  costs  to  “continuation”  configurations. 
In  order  to  do  this,  construct  a  dual  edge  pattern  by  placing  a  point  in  each  cell  of  the  original 
pixel  array,  and  joining  points  if  their  corresponding  pixels  have  some  boundary  in  common. 
The  vertices  of  the  dual  array  can,  in  principle,  be  placed  anywhere  in  their  corresponding  pixels, 
but  in  practice  they  will  have  a  natural  position.  For  example  if  the  pixels  are  constructed  as 
the  Voronoi  polygons  of  a  point  process  then  the  points  of  the  process  will  themselves  be  the 
vertices  of  the  dual  array. 

Our  assumption  that  exactly  three  pixels  meet  at  each  vertex  of  the  original  tessellation 
implies  that  the  dual  edge  pattern  will  be  a  triangulation  of  the  plane.  In  the  case  of  the  square 
pixel  array  the  cost  of  “continuation”  configurations  was  determined  by  considering  a  pattern 
with  a  single  long  straight  edge,  suitably  discretized  to  fit  the  pixel  pattern.  In  the  more  general 
case,  it  is  no  longer  quite  so  clear  how  this  discretization  should  be  performed.  One  natural 
way  to  proceed  is  to  prescribe  that  an  edge  segment  will  be  present  in  the  edge  process  if  and 
only  if  the  corresponding  dual  edge  is  intersected  by  the  straight  line  boundary.  We  assume,  if 


derive  possible  ways  of  penalizing  for  the  “continuation”  configuration  be  given  by  the  presence 
of  the  edges  dual  to  b  and  c  and  the  absence  of  the  edge  dual  to  a.  These  costs  will  be  based 
on  the  general  idea  that  boundaries  should  cost  an  amount  0  per  unit  length;  for  notational 
simplicity  we  shall  assume  henceforth  that  0  =  1,  and  note  that  the  costs  obtained  should  be 
multiplied  by  0  in  the  general  case. 

Let  l  be  a  random  line  in  the  plane,  random  in  a  sense  that  will  be  made  precise  below. 
Let  t t  be  the  length  of  the  intersection  of  £  with  the  triangle  T.  Then  our  first  possible  cost 
for  the  continuation  configuration  be  is 

V\  =  E(£t  |  £  intersects  b  and  c). 

The  motivation  for  this  definition  is  clear.  Summing  It  over  all  triangles  T  gives  the  length  of 
the  line  i ,  neglecting  end  effects,  and,  when  deciding  how  much  to  charge  for  the  configuration 
be,  we  can  only  take  note  of  information  given  by  the  current  clique.  Hence,  by  standard 
statistical  theory,  the  natural  estimator  of  It  is  the  posterior  expectation  Vi  of  It  given  all  the 
information  available. 

A  second  possible  cost  is  given  in  a  slightly  less  transparent  way.  Let  la  be  the  projected 
length  of  the  side  a  on  the  line  £\  this  length  is  to  be  counted  as  negative  if,  as  in  Figure  4.4,  the 
half  of  t  that  intersects  b  and  c  makes  an  angle  of  more  than  |  with  a.  The  second  proposed 
cost  is 


V2  =  E  (\£a  I  t  intersects  b  and  c)  . 


For  such.  9 ,  l  will  intersect  c  and  b  if  and  only  if  it  intersects  c.  The  set  of  lines  at  orientation 
9  that  intersect  c  make  up  a  strip  of  width  c  sin(i?  -  9)  and  so  we  have 


f(9)  oc  csin(f?  -  9)  for  0  <  9  <  B. 
For  —  C  <  9  <  0,  a  similar  argument  yields 

f(9)  a  6sin(C  +  9)  for  0  <  —9  <  C. 


To  calculate  the  constant  of  proportionality,  we  note  that 


/  csin(f?  -  0)d0  +  f  bsin(C  +  9)dO 
Jo  J-c 


f  f 

=  c  I  sin  <f>d<t>  +  b  sin  <f>  d<j> 

Jo  Jo 

=  b  +  c  —  c  cos  B  —  b  cos  C 


=  6  +  c  -  a, 

and  hence  we  have 

{c sin (B  -  9)/(b  +  c  —  a)  0  <  9  <  B 
6sin(C  +  fl)/(6  +  c-a)  -C  <  9  <  0 
0  otherwise. 

To  calculate  Vi,  consider  first  9  >  0.  Given  that  0  =  9  and"lhat  l  intersects  c  and  b,  the 
expected  value  of  is  half  its  value  when  0  =  9  and  l  passes  through  B.  This  length  is,  by 


v  *  u  w «  ■  l  ■  u **  v m  i 


r 


expressed.  Given  that  0  =  0,  we  have  f,  =  a  cos  0,  and  hence 


rB 

V2  =  \a  /  cos  0  /(0)d0 

7-C 

=  (6  +  c  —  a)-1  |y  ^ac cos 0 sin(5  -  0)d0  +  J  \ab  cos&' sin(C  —  0')g?0'|  . 


(4.3) 


The  first  integral  in  (4.3)  is  equal  to 


rB 

/  ^ac{sin  B  +  sin(f?  -  20)}d0  =  ^ac  B  sin  B  =  \&B 
Jo 


where  A  is  the  area  of  the  triangle  T,  and  hence 


V2  =  i(B  +  C)A/(b  +  c-a) 


(4.4) 


Thus  it  is  clear  that  the  formula  for  V2  is  very  simple  and  more  appealing  than  that  for  Vx. 


5.  Regular  Arrays  Revisited 

In  the  last  section  we  defined  two  different  ways  of  obtaining  penalties  for  continuation 
configurations.  One  of  these  was  based  on  the  length  of  the  intersection  of  a  region  in  the  dual 
triangulation  with  a  random  line,  and  the  other  on  the  length  of  projection  of  such  a  region  on 
a  random  line.  It  turned  out  that  the  projection  penalty  gave  a  much  more  elegant  result.  In 
this  section,  we  shall  apply  the  intersection  and  projection  ideas  to  the  regular  square  lattice 
considered  earlier. 

Our  aim  is  to  obtain  costs  for  the  “turn”  and  “continuation”  configurations  as  illustrated 
in  Figure  3.1.  The  dual  of  the  square  lattice  is  itself  a  square  lattice,  and  the  part  of  the  dual 
corresponding  to  a  clique  is  a  single  square  of  side  h  as  in  Figure  5.1. 


will  satisfy 


W)  =  (2  -  V2)-1  sin  (f  -  |tf|) ,  -f  <  9  <  f 

using  simple  calculus  to  find  the  constant  of  proportionality.  The  intersection  length  s  is  equal 
to  h  sec  6 ,  and  hence  the  expected  intersection  length  is 

[  h  sec  9f\{d)d0  =  /  (l  —  |\/2)  1  sec  9  sin  (|  —  8)  d6 
J-i r/4  JO 

/4 

(1  —  tan  9)dd  =  (V2  —  1)_I[0  —  log  sec  d]*JA 
=  (v/2-l)"1  (f  —  2  log  2)  ■ 

Thus  the  “intersection”  penalty  for  a  configuration  of  type  3  in  Figure  3.1  would  be 

(S-Jlog2)-: 
h/(V2-  1)  a  1.06  h. 

To  find  the  “projection”  penalty  for  such  a  configuration,  note  that  the  appropriate  gen¬ 
eralization  of  the  projection  argument  given  in  Section  4  is  to  take  as  penalty  \  (projection  of 
AB  and  DC)  because  both  AB  and  DC  will  be  edges  of  the  irregular  strip  formed  by  the  union 
of  those  dual  pixels  intersected  by  t.  Both  AB  and  DC  have  projection  length  h  cos  6  on  l .  and 
so  the  “projection”  penalty  for  a  configuration  of  type  3  will  be 

/■ir/4  rir/4 

/  hcosdfi(0)d9  =  (l  -  |\/2)  h  /  cosflsin  —  6)  d6 
J-i r/4  JO 

ri r/4 

=  (2  -  V2)~'h  /  {sin  f  -  sin  (2 9  -$)}d9 

Jo 

=  \irh/(V2  -  1)  =  kh  a  0.95  h 
where  k  is  defined  exactly  as  in  Section  3. 


=  (V5-1)-1  f 

Jo 


To  find  the  penalties  for  “turn”  configurations,  the  work  of  Section  4  can  be  used  almost 
directly,  by  noticing  that  both  the  “intersection”  and  “projection”  penalties  will  be  the  same 
as  those  obtained  there,  for  the  case  of  a  line  crossing  the  two  short  sides  of  an  isosceles  right- 
angled  triangle.  Thus  we  set  a  =  h\/ 2,  b  =  c  =  h,  B  =  C  =  ^  and  .4  =  |  in  the  formulas  (4.2) 
and  (4.4). 

We  obtain  as  the  intersection  penalty  for  the  turn  configuration  \h  log  2/(2  -  \/2)  « 
0.59  h  and  for  the  projection  penalty  =  2 -1/2fch  «  0.67  h.  It  is  noteworthy  that 

21 


give  as  the  cost  of  a  “continuation”  as  shown  in  Figure  6.1  the  quantity 

0  /  ffio 

hi  /  cos  0  sin(#o  —  6)d8  j  I  sin(0  -  &o)d9 

rOo  /  r$  o 

=  \h\  J  {sin^o  +  sin(^o  -  2 9)}d6  j  J  sin  9'  d6' 


=  |/i150sin^o/(l  -  cos#o)- 


(6.1) 


The  other  type  of  continuation,  consisting  of  two  edges  of  length  A2,  will  cost  an  amount 
obtained  by  substituting  ho  for  h\  and  t  -  9q  for  9q  in  (6.1),  viz.  \h2{^  -  Oo)  cos#0/(i  —  sin  #o). 

The  general  idea  of  evolving  penalties  for  continuation  configurations  based  on  a  condi¬ 
tional  expected  projection  length  can  of  course  be  extended  to  more  general  polygons  in  the 
dual  tessellation.  The  advantage  of  the  projection  approach  is  that  consistent  penalties  can  be 
written  down  for  cliques  of  different  kinds  that  appear  in  diferent  parts  of  the  same  pattern. 
Thus  the  circular  pixel  grid  shown  in  Figure  4.1 ‘contains  some  vertices  of  degree  3, which  can 
be  dealt  with  using  the  formulas  of  Section  4,  and  some  of  degree  4  whose  dual  polygons  in  the 
dual  tessellation  are,  for  all  practical  purposes,  rectangles  -  which  are  treated  in  this  section. 


